Вздумалось кому-то (не мне) проверять, что матрица между вычислениями не поменялась. Проверять решил просто: считал определитель, сохранял значение и в нужный для проверки момент вычислял определитель опять. Если определитель не изменился, то можно спать спокойно.

На этом математика кончается и начинается песня. Код попадает ко мне — и начинаются глюки на самой свободной ОС, допиленной сумрачным нордическим гением аж до зелёного хамелеона, одиннадцатой версии и второго сервис-пака.

В результате отладки дохожу до такого кода:

                                double a = det(M);
assert(a == det(M));

                              

Ассерт срабатывает. Ладно, добавляю строчку:

                                assert(det(M) == det(M));

                              

Ассерт не срабатывает. Функция всегда возвращает одно и то же значение. Добавляю:

                                double diff = a - det(M);

                              

Результат равен нулю. Причём строго нулю, посмотрел побайтово. Та-ак… Похоже, что имеем вещественное число, в общем случае не равное самому себе. Уже интересно…

                                double a = det(M);
double b = det(M);
assert(a == b);

                              

Ассерт не срабатывает. Пора в дурку…

Ларчик открывался просто. В сопроцессоре все числа обрабатываются в 10-байтовом формате, а double, как известно, 8 байт. Разработчики самого безглючного компилятора возвращали значение в голове стека сопроцессора и забыли нормализовать его до 8 байт. Нормализация происходила только в случае сохранения значения в переменной. Хвост в 2 байта добавлял несколько знаков к мантиссе и вызывал все эти спецэффекты.